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Abstract: We present a theoretical framework fox strain estimation in 
optical coherence elastography (OCE), based on a statistical analysis of 
displacement measurements obtained from a mechanically loaded sample. 
We define strain sensitivity, signal-to-noise ratio and dynamic range, and 
derive estimates of strain using three methods: finite difference, ordinary 
least squares and weighted least squares, the latter implemented for the first 
time in OCE. We compare theoretical predictions with experimental results 
and demonsfrate a ~12 dB improvement in sfrain sensitivity using weighted 
least squares compared to finite difference strain estimation and a ~4 dB 
improvement over ordinary least squares strain estimation. We present 
strain images {i.e., elasto grams) of tissue-mimicking phantoms and excised 
porcine airway, demonstrating in each case clear contrast based on the 
sample's elasticity. 
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1. Introduction 

Disease often causes alteration of the elastic properties of tissue [1,2]. The emerging 
technique of optical coherence elastography (OCE) is aimed at differentiating tissues by 
imaging their microscopic elastic properties [3-6]. In OCE, a tissue sample is subjected to a 
load (stress) and the resulting local displacement is measured using optical coherence 
tomography (OCT) [7]. The spatial derivative of the displacement is the strain [1,2], which is 
indicative of the sample's elastic properties. The strain is displayed in a two- or three- 
dimensional image, known as an elastogram [2]. The limits on resolution and penetration 
depth of displacement measurements in OCE are set by OCT, the underlying imaging 
modality and, thus, lie in the range 1-20 |im and 1-3 mm, respectively. 

Early work in OCE used speckle-tracking algorithms to estimate displacement [3-6]. 
Speckle tracking is complex and, thus, difficult to implement in real time, and has restrictive 
upper and lower limits on displacement, set by speckle decorrelation [4,8-10] and spatial 
sampling in the underlying OCT system, respectively. In phase-sensitive OCE, proposed in 
part to address these issues [11-16], the axial component of displacement is calculated from 
the phase difference between A-scans originating Irom within a single B-scan [11] or Irom 
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successive B-scans [12]. The accuracy and precision of the measured displacement are limited 
by the signal-to-noise ratio (SNR) of the OCT signal [17], and can be corrupted by phase 
wrapping [18]. Since the strain is a differential measure, any errors are amphfied in its 
calculation [19]. To date, there has been no published framework that relates errors in the 
displacement measurement to the calculated strain, making it difficult to optimize OCE, 
assess elastogram image quality and quantitatively compare OCE techniques. 

In this paper, we present a theoretical framework for strain imaging parameters in phase- 
sensitive OCE based on a statistical analysis of the measured displacement. We focus on two 
previously proposed strain estimation methods: finite difference (FD) [3] and ordinary least 
squares (OLS) [20,21]. We define expressions for strain sensitivity, SNR and dynamic range. 
In this analysis, we treat each measurement as a realization of an underlying random variable. 
We compare the performance of both methods using theoretical predictions and experimental 
measurements, and demonstrate the significant sensitivity advantage of OLS sfrain estimation. 
We propose a new modified least-squares sfrain estimation method for OCE, based on 
weighted least squares (WLS). In this method, a weight is assigned to each displacement 
measurement, based on the known OCT SNR at that location. We present elastograms 
generated using each method and demonsfrate an improvement in the sfrain sensitivity of 
WLS strain estimation compared to both FD and OLS strain estimation. Finally, we also 
demonstrate an improvement in strain imaging parameters using Gaussian smoothing. 

The paper is structured as follows: Section 2 describes the measurement of strain with 
OCE and provides a statistical analysis of both FD and OLS strain estimation. In Section 3, 
this framework is used to define strain SNR, sensitivity and dynamic range, and the new WLS 
method for strain estimation is introduced. Section 4 provides details of the experimental 
setup and methods used to validate the theoretical framework. Section 5 presents elastograms 
recorded from phantoms and freshly excised porcine airway. Section 6 presents discussion of 
the results and, finally. Section 7 presents the main conclusions of the work. 

2. Phase-sensitive OCE strain 

Young's modulus, Y, is often used to characterize a material's elastic properties. It is defined 
as the ratio of stress, the force per unit area, to strain, the resulting fractional change in length 
of the sample, and is given by [1] 



where cj is the sfress, is the bulk strain, F is the applied force, A is the surface area, Al is 
the change in length, and /qIs the original length of the sample. The objective of OCE is to 

estimate a sample's depth-resolved elastic properties from OCT-derived displacement data. 
The sfrain can be used as a surrogate for elasticity, if the stress introduced throughout the 
sample is assumed to be uniform [1,2]. However, , defined in Eq. (1), is not depth -resolved; 
it refers to a single measurement of the bulk strain present in a sample. In OCE, the local 
strain, , is a more suitable definition [1-3,20]: 



where Ad is the change in displacement measured over an axial depth range Az , which 
defines the axial resolution of the local strain measurement. Assuming uni-axial loading, if 
Az is the same as the total length of the sample, ^ , then s, and f ^ are equivalent. In general, 
strain measured in elastography refers to local rather than bulk sfrain [1,2]. We will adopt this 
generalization and, therefore, refer to local strain as simply e. 




(1) 
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The goal of strain estimation in OCE is to obtain the strain defined in Eq. (2) fi-om a set of 
measured displacements. Two major factors in determining the strain measurement accuracy 
are the accuracy of the displacement measured from the OCT data and the accuracy of the 
method used to calculate strain irom the measured displacement. 

2.1. Displacement measurement 

A single complex OCT A-scan containing A/^data points may be represented by the set 

/(^•) = {(z,,aP expO«^P))|z, = {i-\)Sz + z, e R,aP e e (3) 

where i is the imaginary unit, z^ are the depths to which the complex OCT measurements 

ap^ exp(z^}^') correspond, dz is the separation between successive z^, and Zj is the first 
depth in the A-scan. In this paper, the superscript j is used to denote a particular A-scan and 
could be, for example, one of many making up a B-scan, or two A-scans acquired in the same 

lateral position in the sample in successive B-scans. The phase, , is in general random, 

although a map of the displacement of scattering centers, occurring between the acquisition of 

A-scans 7*^' and /^'^ can be found, without the loss of generality, from the phase difference 
between complex OCT signals as 



D= (z,,af,) 



Ann Ann 

i i=i 



where X is the mean source wavelength and n is the refractive index of the material. 
Displacements deriving from a phase change of greater than ±71 are ambiguous. This is a 
limitation referred to as phase wrapping [18] and sets the upper limit of measurable phase 
difference, without the need for complicated phase imwrapping algorithms. 

The minimum measurable phase difference is limited by the standard deviation of the 
random variable for which A^, is a realization [17]. Under the assumption that the SNR of the 

detected OCT signal intensity SNRg^^.^ satisfies, SNRq^t. » 1 , it can be shown that the phase 
sensitivity, (t^^. , is given by [17]: 



(5) 



OCT, 



For example, considering a typical SNRq^j. of 20 dB, combining Eqs. (4) and (5) gives a 
displacement sensitivity of 4.7 nm, assuming A = 800 nm. It should be noted that if 
SNRqqj » 1 , the phase difference is unbiased; thus, the expected value of the displacement 

random variable corresponds to the true sample displacement [22]. 
2.2. Strain estimation 

In OCE, two methods have been proposed to estimate strain from the measured displacement 
map; FD strain estimation [3] and OLS strain estimation [20,21]. To determine the accuracy 
and precision of these methods, we define expressions for their mean and variance. 

2.2.1. Finite difference strain estimation 

Strain estimation using the FD method was previously proposed in the related field of 
ultrasound elastography [23], and later applied to OCE [3]. FD strain estimation, sf , is 
defined as 
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^i+m-\ ' 



Az ' 



(6) 



where Az = mSz . As the strain is obtained from the hnear combination of two random 
variables, its variance can be defined as [24]: 



,<+<.-,-2cov(A-,A-..-i) 



(7) 



where E/ , Z), and are random variables of which e{ , d. and fi?/+,„_] are reahzations. 

Assuming the measured displacements are from non-overlapping positions, they can be 
assumed to be statistically independent [24] and the covariance term in Eq. (7) can be set to 
zero. The expression is further simplified assuming the displacement measurements are 
homoscedastic, i.e., the variance is the same for D. and Df^^_^ , and is written simply as cj^ : 



2c^ 

' Az' 



(8) 



From Section 2.1, the displacement estimation is unbiased. Therefore, FD strain estimation is 
also imbiased; the mean of the random variable Ef corresponds to the true underlying strain: 



lip/ 



(9) 



2.2.2. Ordinary least squares strain estimation 

An alternative method to estimate strain is ordinary least squares (OLS). In contrast to the FD 
method, in the OLS method use is made of more than two data points to obtain greater 
accuracy. In ultrasound elastography, the OLS method has been shown to provide more 
accurate strain estimation than the FD method [19]. We aim here to provide such a 
comparison in OCE for the first time. By restricting z to the interval z^ < z < z,. + Az , we can 

rewrite Eq. (2), imphcitly defining £: 



d = s{z — z._^') + c. 



(10) 



Given a set of depth values and associated displacement measurements, d, within the interval 
defined above, the strain is determined by minimizing the sum of squared differences, R, 
between the measured displacements and the displacements given by Eq. (10), i.e., by 
minimizing the following expression: 



i+m-l 



(11) 



A complete derivation of OLS is provided elsewhere [25]. Following this derivation, the strain 
can be estimated using the OLS method according to 



£t = 



i+m-\ 



I t +m- \ ' 

^1 J^{Zj-Z,^,f - Y^{zj-z,_,) 

V j=i A j=> J V J=' 



(12) 
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where m is the number of data points used in the strain calculation. Using this expression, we 
obtain an estimation of strain at each data point, i, within an A-scan. Equation (12) can be 
simphfied by reducing the number of summation operators such that 



'^o(2;-Vl)-'^l 



dj. 



where 



(13) 



(14) 



Substituting Zj - z,_j = (y - i + 1)& in Eqs. (13) and (14) allows s" to be expressed as 



i+m-l/' / . -.1X0 ^ 

«'oO-^ + i)&-«'i 



K:qK2 — Ki 



(15) 



Equation (14) can now be simplified using mathematical induction to remove the summation 
operators: 



i+m-l ^ 



Substituting Eqs. (16)-(18) in Eq. (15) and simphfying, we obtain 



(16) 



(17) 



(18) 



(19) 



To allow expressions for the mean and variance of £° to be derived, we assume that the set of 
displacements has the same variance, i.e., homoscedasticity, and use the scaling and summing 
properties of independent Gaussian random variables. The mean of E° , the random variable 

of which is a realization, can be expressed as 



^E' = Z 



"■-'^ 6(2(7 -z)-OT + l) 



Sz(m —m) 



(20) 



As the random variable Z)^. is unbiased (see Section 2.1), its expected value, //^ , equals the 

true displacement. Replacing z-z^_^ with {j-i + \)5z in Eq. (10), and then substituting this 
into Eq. (20), and simplifying using mathematical induction, we obtain 



6(2(7-0-m + l)V X 

^E' = Z -^\{e{j-i + \)dz + c): 

Sz{m -m) ) 



(21) 
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Equation (21) con&ms that the mean of E° equals the true underlying strain. Similarly, the 
variance of E° may be defined as 



6(2(7-0 -« + l) 



ol = '"7 (22) 

Sz {m —m) Az m 



I2cl 12cj 



Sz(m —m) 
3. Strain imaging parameters 

Having derived the mean and variance of the strain estimated using both FD and OLS 
methods, we now define the strain sensitivity, SNR and dynamic range in terms of these 
parameters. The strain sensitivity is the minimum detectable strain in an elastogram. 
Analogously to the definition of phase difference sensitivity in Eq. (5), the strain sensitivity 
for FD and OLS strain estimation, s^/io , is defined as the standard deviation of the strain 

[26]: 

In ultrasound elastography, the strain SNR, SNR^/io , is defined as the ratio of the mean to the 
standard deviation of the strain [26]. We choose the same definition for strain SNR in OCE: 

SNRf,„=^^. (24) 

Strain dynamic range, DR^^,, , is the ratio of the largest to the smallest measurable strain [26]. 

Maximum strain in OCE occurs when the maximum measureable change in displacement, 
'!^<^max5 occurs over one strain axial resolution element. To avoid the ambiguity caused by 
phase wrapping, the upper limit of displacement considered in this paper corresponds to a 
phase difference of n radians. Minimum strain is defined by ^ //„ , allowing the dynamic 

range to be defined as 




3.1. Strain imaging parameters using the FD method 

Substituting the mean and variance of the strain estimated using the FD method, defined in 
Section 2.2.1, into the definitions of the strain imaging parameters describes the performance 
of FD strain estimation. The FD strain sensitivity, S^f , SNR, SNR^f , and dynamic range, 

DR^f , are foimd to be 

^B.-4^; (26) 

Az 

SNR^,=-jL — ;and (27) 
' V2(To 

DR^f-^T^- (28) 



32f7o« 
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3.2. Strain imaging parameters using the OLS method 



Similarly, substituting the mean and variance of the strain estimated using OLS, defined in 
Section 2.2.2, in the expressions for the strain imaging parameters defined above allows us to 
assess the performance of this method. Using this approach, the OLS strain sensitivity, , 

SNR, SNR^o , dynamic range, DR^o , are defined as 



s^^-f^; (29) 

Az^m 
\x Az^Jm 

SNR „ = ; and (30) 



' D 




(31) 



3.3. Theoretical predictions 



In Fig. 1, we present theoretical predictions obtained using the expressions derived in Sections 
3.1 and 3.2 for FD and OLS strain estimation. The OCT SNR and source mean wavelength 
were assumed to be 20 dB and 835 nm, respectively. The variation in strain sensitivity with 
strain axial resolution, presented in Fig. 1(a), demonstrates the increasing sensitivity 
advantage of OLS over FD strain estimation as the resolution degrades. An improvement in 
sensitivity of ~8 dB is achieved using OLS estimation at the lowest simulated axial resolution 
of Az= 120 ^m. 

In Fig. 1(b), the SNRs of the strain using FD and OLS estimation are compared. The strain 
was chosen to be 1200x10 *■ , corresponding to the maximum change in displacement, Aiimax, 
over the lowest strain axial resolution considered, i.e., 120 |a,m. As discussed in Sections 3.1 
and 3.2, the expected value of strain using both methods is unbiased and, therefore, 
corresponds to the true strain. As the SNR is the ratio of mean strain to strain sensitivity, the 
difference between the FD and OLS methods is determined by the strain sensitivity. This is 
confirmed in Fig. 1(b), where the SNR of the OLS estimate is also improved by ~8 dB for Az 
= 120 |a,m. Figures 1(a) and 1(b) highlight a fiindamental trade-off in strain estimation in OCE, 
namely, improved strain sensitivity and SNR, using both methods, at the expense of axial 
resolution. 

In Fig. 1(c), the SNR of FD and OLS strain estimations versus the change in displacement, 
Afi?, measured in one strain axial resolution element, is presented. The resolution was fixed at 
120 nm, corresponding to the lowest value shown in Figs. 1(a) and 1(b). For each Ad, the 




0 20 40 60 80 100 120 0 20 40 60 80 100 120 0 50 100 150 



Strain resolution, \xm Strain resolution, nm Ad, nm 



Fig. 1 . Theoretical predictions of strain estimation using FD (dashed red) and OLS (solid blue) 
methods: (a) sensitivity; and (b) SNR versus strain axial resolution, (c) SNR versus change in 
displacement, Ao*, in one resolution element. 
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SNR of the OLS estimate is higher than that of the FD estimate. Figure 1(c) highlights the 
importance of maximising Ad in each strain resolution element. 

3.4. Weighted least squares strain estimation 

According to the Gauss-Markov theorem [25], OLS provides the best linear unbiased strain 
estimation if the errors have the same variance {i.e., assuming homoscedasticity). However, in 
phase-sensitive OCE, the variance of each displacement measurement is dependent on the 
OCT SNR (Eq. (5)). To address this and to improve the strain estimation, we adopt a weighted 
least-squares (WLS) method and assign a weight to each displacement measurement based on 
the known variance of that measurement. A strain estimation method based on WLS has been 
proposed in ultrasound elastography [27], however, we report its first use in OCE. The 

derivation of WLS strain, , follows the same procedure as for £° , enabling the following 

expression to be defined: 



£ = 



^i+m-l ^ f i+in-\ ^ f \ f i+m-l \ 

Z 2] w.(z.-z._,)(/. - ^w/z.-z,_i) J^Wjdj 



]f 



fi+m-l 



(32) 



■) 



where the weighting variable, Wj , is chosen to be 



(33) 



In Eq. (33), is determined by combining Eqs. (4) and (5). Analogously to £° , the strain 

axial resolution, Az , and the change in displacement over this range. Ad, are determined by 
the number of data points, m, over which is calculated. As this method incorporates more 
information about the measured displacement, it has the potential to provide more accurate 
strain estimations than either the FD or OLS methods. Note that because w. depends on 

SNRg(.j. , which is non-uniform and not known a priori, we are unable to derive for WLS 

strain estimation simple analytical expressions for strain imaging parameters from Eq. (32). 
Instead, we demonstrate its superiority experknentally in Section 5. 

4. Experimental setup and procedure 

A schematic diagram of the SD-OCE system is presented in Fig. 2(a). The system employs a 
superluminescent diode source with mean wavelength of 835 nm and 55.4 nm bandwidth, 
illuminating the sample with 10 mW of optical power. The theoretical lateral resolution is 15 
(jm and the sensitivity was measured to be 108 dB for an exposure time of 98 |is. The 
spectrometer consists of a diffraction grating with 1800 hnes/mm, a lens with a focal length of 
50 mm and a line-scan camera, consisting of 4096 pixels. The sensitivity roll-off was 
measured to be 4 dB/mm. The A-scan acquisition time was 100 (is. Lateral scanning over a 4 
mm range was performed using a galvanometer mirror, driven by a sawtooth waveform. For 
each elastogram, 1000 A-scans were recorded in an acquisition time of 0.1 s. 

Samples were mechanically loaded using a piezoelectric ring actuator described 
previously [20] that performed imaging and compression from the same side of the sample. 
Samples were preloaded by compression against a rigid upper boundary to ensure full contact 
with the sample surface. A 5 Hz square wave synchronized with the lateral scanning was 
apphed to the actuator, as illustrated in Fig. 2(b), to produce a displacement amplitude of 0.15 
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|im, corresponding to the upper limit set by phase wrapping. The resulting sample 
displacement was calculated from the phase difference between A-scans at the same lateral 
position acquired in successive B-scans. We note that by using this procedure the phase 
sensitivity is not affected by the lateral oversampling in the x-direction, as is the case when the 
phase difference is measured between A-scans acquired within one B-scan [17]. However, 
when performing 3-D OCE measurements, scanning in the y-direction degrades phase 
sensitivity, and oversampling in this axis is required to minimize this reduction in phase 
sensitivity. 

In Fig. 2, the phase difference, measured using a stationary glass slide in the sample arm, 
is presented for each of 100 B-scans, both without. Fig. 2(c), and with. Fig. 2(d), lateral 
scanning. In both cases, the OCT SNR was measured to be 50 dB. In the non-scanning case, 
(T^p was measured to be 4 mrad, within 1 mrad of the value given by Eq. (5), and 

corresponding to a displacement sensitivity of 0.2 nm. In the scanning case, cr^^ was 

measured to be 25 mrad, corresponding to a displacement sensitivity of 1.2 nm. This is similar 
to the phase sensitivity reported previously in OCT systems with galvanometer-based lateral 
scanning [28]. The increased phase noise is due to galvanometer jitter between B-scans and is 
not related to oversampling, as no scanning is performed in the j-direction. 

For the elastograms presented in Section 5.2, the variance of the measured displacement 
was reduced by averaging 100 B-scans at each square wave half-cycle prior to calculation of 
the phase difference. The averaging was performed in the complex plane, before extracting 
the phase, in order to avoid the systematic error caused by the 2n phase modulus [18,29]. 

Both phase and strain sensitivity are determined by the OCT SNR (Eq. (5)), which is 
strongly affected by speckle and signal attenuation. To improve the strain sensitivity, 
Gaussian smoothing (GS) was applied by convolving the elastogram with a Gaussian kernel 
with frill- width at half maximum (FWHM) of 20 nm x 20 ^m. 
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Fig. 2. (a) Schematic diagi'am of the OCE system; RB: rigid boundary; RA: ring actuator; SM: 
scamiing min'or; BS: beam splitter; RM: reference mirror; SLD: superluminescent diode; (b) 
Saw-tooth pattern of lateral scanning beam synchronized with motion of ring actuator, (c) The 
measured phase difference with no lateral scanning and (d) with lateral scamiing. 

4.1. Phantom fabrication 

Tissue-mimicking phantoms with approximate dimensions (length x width x height) of 10 
mm X 10 mm x 2 mm were fabricated using room-temperature vulcanizing (RTV) silicone as 
the bulk medium and titanium dioxide (TiOi) as optical scatterers. Three phantom designs 
were fabricated. Phantom 1 was designed to be optically and mechanically homogeneous and 
relatively soft. It was fabricated from sofr silicone (Elastosil® P7676, Wacker, Germany) and 
contained 1.5 mg/ml of Ti02 particles. Its elastic modulus was measured using an Instron 
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materials testing system to be 25 kPa. Phantom 2 was designed to be a soft bulk medium 
containing hard inclusions. It was fabricated using the same silicone used in Phantom 1 . The 
inclusions were fabricated using a hard silicone (Elastosil® RT 601, Wacker, Germany), with 
measured elastic modulus of 4 MPa. Inclusions of measured average diameter of 0.4 mm were 
cut from a block by hand using a scalpel under a dissecting microscope. The inclusions 
contained 3.5 mg/ml Ti02 particles. They were embedded in the soft surrounding silicone 
using a two-stage curing process. A 1.5 mm-thick layer of soft silicone was cured in a mold. 
The inclusions were placed on the surface of the soft layer using a hypodermic needle. A 
second layer of soft silicone, from the same batch as the first layer, was poured on top of the 
inclusions and allowed to cure, encasing them in the soft medium of total thickness 2 mm. 
Phantom 3 was designed to be a bilayer phantom with hard material of thickness 0.6 mm on 
top of soft material of thickness 1.4 mm. The top layer was fabricated using RT 601 silicone 
and contained 1.5 mg/ml Ti02 scatterers. The bottom layer was fabricated using P7676 
sihcone and contained 3.5 mg/ml Ti02. The elastic moduli of the hard and soft layers were 
measured to be 4 MPa and 25 kPa, respectively. 

5. Experimental results 

5. 1. Comparison between theoretical predictions and experimen t 

In this section, experimental results are presented and compared to theory. Experiments were 
performed using Phantom 1 with uniform elastic properties so that the strain introduced was 
constant for all positions in the sample. To maximize phase sensitivity, no lateral scanning 
was performed. The theoretical and measured sfrain sensitivity, mean and SNR are presented 
in Fig. 3. The measured parameters were calculated from experimental data for FD (red dots), 
OLS (blue dots), WLS (green dots) and Gaussian smoothed weighted least squares (GS-WLS, 
black dots) strain estimation methods. For each method, the mean and standard deviation of 
100 calculations from consecutive axial positions within the sample were obtained, for strain 
resolutions from zero to 120 |j,m. The strain sensitivity and SNR were calculated from the 
experimental data using Eqs. (23) and (26). The theoretical predictions of the strain sensitivity 
and SNR are shown for the FD (solid red) and OLS (solid blue) methods. To best compare 
experiment and theory, the standard deviation of the experimentally measured displacement 

was substituted for cjj^ in the expressions for 'E{ and E° defined in Sections 3.1 and 3.2. In 

Fig. 3, close agreement is seen between experiment and theory for both FD and OLS strain 
estimation methods. As expected, for both strain sensitivity (Fig. 3(a)) and strain SNR (Fig. 
3(c)), the WLS method provides better performance than FD (by ~12 dB) and OLS estimates 
(by ~4 dB). The results obtained using GS-WLS strain estimation are also presented in Figs. 
3(a)-3(c). The smoothing results in a ~3 dB improvement in sfrain sensitivity and SNR in 
comparison to WLS estimation. The mean strain calculated for each method is presented in 




0 20 40 60 80 100 120 0 20 40 60 80 100 120 0 20 40 60 80 100 120 



Strain resolution, |.im Strain resolution, Strain resolution, |.im 

Fig. 3. (a) Strain sensitivity; (b) mean strain; and (c) strain SNR of Phantom 1. Experimental 
results are presented for FD (red dots), OLS (blue dots), WLS (green dots), and GS-WLS 
(black dots) strain estimation. Theoretical results are presented for FD (red line) and OLS (blue 
line) strain estimation. 



#168336 - $15.00 USD Received 10 May 2012; revised 28 Jun 2012; accepted 12 Jul 2012; published 17 Jul 2012 
(C) 2012 OSA 1 August 2012 / Vol. 3, No. 8 / BIOMEDICAL OPTICS EXPRESS 1875 



Fig. 3(b). For strain resolutions >40 \im, the mean strain converges to -38 dB for all methods. 
This value was, therefore, used in the simulations. Below -40 (xm, the mean strain measured 
with all methods varies by ~3 dB, suggesting that below this threshold the measured strain is 
unrehable. 

5.2. Elastograms 

In this section, elastograms and the corresponding structural OCT images are presented for 
Phantoms 2 and 3, and for an excised porcine airway sample. For each image, the scale bars 
correspond to physical depth, assuming a refractive index of 1 .4, the strain axial resolution is 
75 nm, and the lateral resolution is 15 (xm. Imaging and compression were introduced from 
the top in all elastograms. Each sample was tilted with respect to the beam axis to minimize 
spectral reflection from the air-sample interface. To maintain uniform compression, the ring 
actuator compression plate and the upper rigid boundary were tilted by the same angle in both 
X and y planes. In each elastogram, the strain SNR is plotted. Each pixel in an elastogram 
corresponds to the local strain divided by the noise floor, which was determined from the 
standard deviation calculated in an area (50 (xm x 50 (xm) with uniform intensity in the OCT 
structural image. 

5.2.1. Phantom 2: Soft with hard inclusions 

In Fig. 4(a), the OCT structural image of a central portion of Phantom 2 is presented. A 
rectangular-shaped inclusion is visible commencing at -200 |xm below the surface. In Figures. 
4(b)-4(e), the corresponding elastograms are presented for (b) FD, (c) OLS, (d) WLS and (e) 
GS-WLS strain estimation methods. Using FD strain estimation (Fig. 4(b)), it is difficult to 
distinguish the inclusion due to the low strain SNR. The result obtained with OLS, Fig. 4(c), 
is significantly improved, as expected, and the inclusion is clearly distinguishable. The strain 
SNR is further improved using WLS, as evident in Fig. 4(d). 

Large striations or bands are evident in the elastograms presented in Figs. 4(b)-4(d). The 
origin of this banding is discussed in Section 6. Its effect can be reduced using Gaussian 
smoothing, shown in Fig. 4(e). The tradeoff inherent in Gaussian smoothing is the 
improvement in SNR at the expense of spatial resolution. In Fig. 4(e), the axial resolution is 
marginally reduced from 75 (xm to 78 (xm, but the lateral resolution is degraded from 15 (xm to 
25 |xm. Figures 4(f)-4(j) further illustrate the performance of the different strain estimation 
methods, showing the variation of OCT and strain SNR over the same lateral range as the 
images, at the depth indicated by the blue arrow in Fig. 4(a). The high spatial frequency noise 
in FD, OLS and WLS strain estimation is visible in Figs. 4(g)-4(i), as is its reduction in Fig. 
4(j)- 

The low OCT SNR in Fig. 4(a) below the hard inclusion is a shadow artifact caused by 
signal attenuation in the highly scattering inclusion. High strain is observed in this location in 
all elastograms. From a mechanical standpoint, this is expected due to the increased stress 
caused by the inclusion. However, as the OCT SNR decreases, the probability distribution of 
the measured phase difference tends toward uniform [22]. This reduction in the mean phase 
difference causes the measured strain to increase below the inclusion. This is an artifact in 
strain estimation and is visible in all elastograms in Fig. 4. It could potentially be removed 
using additional thresholding to set a lower limit on OCT SNR at which strain can be 
measured rehably. 
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Fig. 4. Phantom 2: (a) OCT structural image; and (b) FD; (c) OLS; (d) WLS; and (e) GS-WLS 
elastograms. In (f)-(j), corresponding lateral traces are shown for the depth indicated by the 
blue arrow in (a). 

5.2.2. Phantom 3: Hard/soft bilayer phantom 

hi Fig. 5(a), the OCT structural image of Phantom 3 shows the layer boundary at a depth of 
-600 i^m. In Fig. 5(b), the corresponding elastogram, generated using GS-WLS, shows that 
the layers can also be distinguished based on strain. The low strain SNR in the hard first layer 
suggests that it is displaced predominantly as a bulk, with most of the stress being dissipated 
in the underlying soft layer, as expected [20,30]. The banding effect described in Section 5.2.1 
is also visible in Fig. 5(b), and will be discussed in Section 6. It can be seen in Fig. 5(b) that 
the strain in the second layer is non-uniform. This is due to two factors. Firstly, the silicone, 
like many soft tissues, is incompressible and therefore has a Poisson's ratio of -0.5. In this 
case, axial deformation detected using phase-sensitive OCE, is accompanied by a 
corresponding lateral deformation. In the softer second layer, near the boundary with the 
stiffer first layer, the lateral deformation of the sample is restricted by adhesion to the stiffer 
layer. This causes lower strain in the softer layer near the boundary, resulting in a higher 
apparent Young's modulus. Secondly, the lower axial resolution in the elastogram reduces the 
sharpness of interfaces and contributes to the strain gradient in the second layer in Fig. 5(b). 
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5.2.3. Porcine airway 

To demonstrate an elastogram of tissue, we performed measurements on fresh ex vivo porcine 
airway harvested from a heahhy animal. The airway was sectioned into a flat 10 mm x 10 mm 
X 2 mm sample. The sample was positioned such that the inner wall of the airway was in 
contact with the actuator compression plate. An OCT image of the airway is presented in Fig. 
5(c). Three distinct layers are visible and readily attributable to the expected morphology. The 
highly scattering, most superficial layer is the mucosa; the middle low- scattering layer is the 
cartilage; and the deepest layer is the adventitia. The layers visible in Fig. 5(c) correspond 
well with those reported in previous OCT airway studies [31,32]. In Fig. 5(d), the GS-WLS 
elastogram is presented (generated in the same manner as for Figs. 4(e) and 5(b)). The three 
layers are apparent in the elastogram. Both the innermost mucosal layer and the outermost 
adventitial layer have relatively high strain, in comparison to the stiff cartilage, which has 
relatively low strain. The high stiffness of the cartilage is expected due to its composition, 
largely of collagen fibers [33]. 




Fig. 5. (a) OCT structural image; and (b) GS-WLS elastogram of Phantom 3. (c) OCT 
structural image; and (d) GS-WLS elastogram of a thin section of excised porcine airway with 
layers as labeled in (c). 



6. Discussion 

In any imaging technique, a framework characterizing performance, defining parameters such 
as spatial resolution, sensitivity, SNR, and dynamic range, enables assessment and 
optimization of image quality. Such a framework has played an important role in the closely 
related technique of ultrasound elastogi^aphy [19,24,26]. In this paper, we have presented the 
first such framework for OCE in considering the performance of three strain estimation 
methods: two previously proposed, based on FD and OLS strain estimation, and a method that 
we have proposed for the first time in OCE, WLS strain estimation. Importantly, the 
application of these methods in OCE is distinct to their application in ultrasound elastography 
because of the different noise characteristics and processing steps used in each case. For 
example, in ultrasound elastography, the variance of the measured displacements is 
determined by parameters such as the transmitter bandwidth and the correlation coefficient 
[26], whereas in phase-sensitive OCE, statistical optics determines this variance [22]. 

In the elastograms presented, a banding effect is visible, i.e., the modulation between high 
and low strain in adjacent axial regions of elastograms. We demonstrated that Gaussian 
smoothing reduces this effect, at the expense of degradation of the spatial resolution. Future 
directions will include the investigation of more sophisticated filters to mitigate this 
degradation, e.g., edge-preserving bilateral filters [34]. These banding artifacts are similar in 
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appearance to those reported in ultrasound elastography [2]. Their origin has been attributed 
to cyclic bias error in displacement measurements and correlated noise patterns arising from 
the use of oversampling. This is the first time such artifacts have been reported in OCE, and 
further work is required to confirm their origin and develop methods for their mitigation. We 
anticipate that these artifacts are related to the change in the probability distribution of the 
phase difference measurements when the condition SNRg^^j. » 1 is not satisfied. Under these 
conditions, Eq. (5) is not valid and the phase difference distribution is uniform rather than 
Gaussian [22]. This introduces a systematic bias error in the measured displacement, as the 
mean displacement, //^ , is zero rather than the true displacement. This bias influences the 

estimated strain and may explain the observed banding effect. Low SNR largely results from 
attenuation as light propagates through a sample, which could explain why the banding effect 
is more prominent at greater depths in the elastograms presented in Figs. 4 and 5. 

7. Conclusions 

We have presented a theoretical framework for strain estimation in OCE based on FD, OLS 
and WLS methods. We have defined, for the first time, expressions for strain sensitivity, SNR 
and dynamic range. For FD and OLS strain estimation, we demonstrated both close agreement 
between theory and experiment, and the advantage of OLS strain estimation. Our 
experimental results indicate that strain estimation using WLS provides more accurate results 
than using FD (by ~12 dB) or OLS (by ~4 dB). We demonstrated that Gaussian smoothing 
conveyed an additional improvement in strain sensitivity at the expense of lateral resolution. 
The framework presented provides a quantitative means to relate phase sensitivity to 
elastogram image quality in phase-sensitive OCE. We expect that this framework will allow 
quantitative comparisons to be made between different strain estimation methods, as well as 
different OCE experimental procedures, greatly facilitating the future development of OCE. 
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